
from numpy import *
from cmath import exp
from plotter import Plotter

from scipy.spatial import KDTree

def Init_Hops_NaxCoO2_Frank_L3():
  """
  Initializes hopping parameters for NaxCoO2 with 3 orbitals,
  a1g, eg, eg, included.
  """

  thop = zeros((3,3,9), float)

#
# Notations:
#   1. Plain index notation (like in the notes).
#   2. Absolute normalized coordinates.
#   3. Direct coordinates in vectors: 
#        a1 = {1/2, -sqrt(3)/2}, 
#        a2 = {1/2,  sqrt(3)/2},
#        a3 = { 0, 0, 1}
#
#      (Frank's data).
#
# t1 = t_{ 1/2, sqrt(3)/2 } = t_{0, 1, 0}
  thop[:,:,0] = \
   [[ 137.0,    45.1,    57.0], \
    [  27.2,   -75.7,    33.9], \
    [  67.4,    22.4,   -43.4]]

# t2 = t_{ 1, 0 } = t_{1, 1, 0}  
  thop[:,:,1] =  \
   [[ 137.2,   -71.7,   -10.3 ], \
    [ -71.7,   -27.4,    -5.8 ], \
    [  10.3,     5.8,   -92.0 ] ]

# t3 = t_{ 1/2, -sqrt(3)/2 } = t_{1, 0, 0}
  thop[:,:,2] =  \
   [[ 137.0,    27.2,   -67.4 ], \
    [  45.1,   -75.7,   -22.4 ], \
    [ -57.0,   -33.9,   -43.4 ] ]
 
# t4 = t_{ 3/2, sqrt(3)/2 } = t_{ 1, 2, 0 }
  thop[:,:,3] =  \
   [[ -19.5,    -6.0,    10.4 ], \
    [  -5.8,     8.5,    12.4 ], \
    [   9.9,    12.4,    -5.8 ] ]

# t5 = t_{ 0, sqrt(3) } = t_{ -1, 1, 0 }
  thop[:,:,4] =  \
   [[ -19.5,    11.5,     0.0 ], \
    [  12.0,   -12.9,     0.0 ], \
    [   0.0,     0.0,    15.6 ]] 

# t6 = t_{ 3/2, -sqrt(3)/2 } = t_{ 2, 1, 0 }
  thop[:,:,5] =  \
   [[ -19.5,    -5.8,    -9.9 ], \
    [  -6.0,     8.5,   -12.4 ], \
    [ -10.4,   -12.4,    -5.8 ]]

# t7 = t_{ 1, sqrt(3) } = t_{ 0, 2, 0 }
  thop[:,:,6] =  \
   [[ -22.3,    18.2,    31.7 ], \
    [  18.3,   -15.6,   -31.4 ], \
    [  31.7,   -29.6,   -50.9 ]]

# t8 = t_{ 2, 0 } = t_{ 2, 2, 0 }
  thop[:,:,7] =  \
   [[ -22.2,   -36.5,     0.0 ], \
    [ -36.5,   -68.5,     0.9 ], \
    [   0.0,    -0.9,     2.0 ]]

# t9 = t_{ 1, -sqrt(3) } = t_{ 2, 0, 0 }
  thop[:,:,8] =  \
   [[ -22.3,    18.3,   -31.7 ], \
    [  18.2,   -15.6,    29.6 ], \
    [ -31.7,    31.4,   -50.9 ]]

# Proper sign for hopping parameters
  thop = -thop / 1000

# *************** DEBUG ************
#  thop[:,:,3:] = 0.0
#  thop[1:,1:,0:3] = 0.0
#  thop[1:,0,0:3] = 0.0
#  thop[0,1:,0:3] = 0.0
#  thop[2:,0:2,0:3] = 0.0
#  thop[0:2,2:,0:3] = 0.0

  eonsite = array( [0.0, -36.2, -36.2] ) / 1000
#  eonsite = array( [0.0, -150.0, -150.0] ) / 1000

  return thop, eonsite

def Init_Hops_L3_deb():
  """
  Initializes hopping parameters for NaxCoO2 with 3 orbitals,
  a1g, eg, eg, included.
  Playground for parameters.
  """

  thop = zeros((3,3,9), float)

#
# Notations:
#   1. Plain index notation (like in the notes).
#   2. Absolute normalized coordinates.
#   3. Direct coordinates in vectors: 
#        a1 = {1/2, -sqrt(3)/2}, 
#        a2 = {1/2,  sqrt(3)/2},
#        a3 = { 0, 0, 1}
#
#      (Frank's data).
#
  t1 = 177.0
  t2 = -12.5
  t3 = -14.3
# t1 = t_{ 1/2, sqrt(3)/2 } = t_{0, 1, 0}
  thop[:,:,0] = \
   [[    t1,    45.1,    57.0], \
    [  27.2,   -75.7,    33.9], \
    [  67.4,    22.4,   -43.4]]

# t2 = t_{ 1, 0 } = t_{1, 1, 0}  
  thop[:,:,1] =  \
   [[    t1,   -71.7,   -10.3 ], \
    [ -71.7,   -27.4,    -5.8 ], \
    [  10.3,     5.8,   -92.0 ] ]

# t3 = t_{ 1/2, -sqrt(3)/2 } = t_{1, 0, 0}
  thop[:,:,2] =  \
   [[    t1,    27.2,   -67.4 ], \
    [  45.1,   -75.7,   -22.4 ], \
    [ -57.0,   -33.9,   -43.4 ] ]
 
# t4 = t_{ 3/2, sqrt(3)/2 } = t_{ 1, 2, 0 }
  thop[:,:,3] =  \
   [[    t2,    -6.0,    10.4 ], \
    [  -5.8,     8.5,    12.4 ], \
    [   9.9,    12.4,    -5.8 ] ]

# t5 = t_{ 0, sqrt(3) } = t_{ -1, 1, 0 }
  thop[:,:,4] =  \
   [[    t2,    11.5,     0.0 ], \
    [  12.0,   -12.9,     0.0 ], \
    [   0.0,     0.0,    15.6 ]] 

# t6 = t_{ 3/2, -sqrt(3)/2 } = t_{ 2, 1, 0 }
  thop[:,:,5] =  \
   [[    t2,    -5.8,    -9.9 ], \
    [  -6.0,     8.5,   -12.4 ], \
    [ -10.4,   -12.4,    -5.8 ]]

# t7 = t_{ 1, sqrt(3) } = t_{ 0, 2, 0 }
  thop[:,:,6] =  \
   [[    t3,    18.2,    31.7 ], \
    [  18.3,   -15.6,   -31.4 ], \
    [  31.7,   -29.6,   -50.9 ]]

# t8 = t_{ 2, 0 } = t_{ 2, 2, 0 }
  thop[:,:,7] =  \
   [[    t3,   -36.5,     0.0 ], \
    [ -36.5,   -68.5,     0.9 ], \
    [   0.0,    -0.9,     2.0 ]]

# t9 = t_{ 1, -sqrt(3) } = t_{ 2, 0, 0 }
  thop[:,:,8] =  \
   [[    t3,    18.3,   -31.7 ], \
    [  18.2,   -15.6,    29.6 ], \
    [ -31.7,    31.4,   -50.9 ]]

# Proper sign for hopping parameters
  thop = -thop / 1000

# *************** DEBUG ************
#  thop[:,:,3:] = 0.0
#  thop[1:,1:,0:3] = 0.0
#  thop[1:,0,0:3] = 0.0
#  thop[0,1:,0:3] = 0.0
#  thop[2:,0:2,0:3] = 0.0
#  thop[0:2,2:,0:3] = 0.0

  eonsite = array( [0.0, -36.2, -36.2] ) / 1000

  return thop, eonsite


def Init_Hops_NaxCoO2_Korshunov_L3_x33():
  """
  Initializes hopping parameters for NaxCoO2 with 3 orbitals,
  a1g, eg, eg, included.
  Parameters are from Korshunov et al.  (x = 0.33)
  """

  thop = zeros((3,3,9), float)

#
# Notations:
#   1. Plain index notation (like in the notes).
#   2. Absolute normalized coordinates.
#
#      (The data is from Korshunov et al.).
#
# t1 = t_{ 1/2, sqrt(3)/2 }
  thop[:,:,0] = \
   [[ 123.0,   -44.0,   -77.0], \
    [ -44.0,   -69.0,    37.0], \
    [ -77.0,    37.0,   -26.0]]

# t2 = t_{ 1, 0 }
  thop[:,:,1] =  \
   [[ 123.0,    89.0,     0.0 ], \
    [  89.0,    -5.0,     0.0 ], \
    [   0.0,     0.0,   -90.0 ] ]

# t3 = t_{ 1/2, -sqrt(3)/2 }
  thop[:,:,2] =  \
   [[ 123.0,   -44.0,    77.0], \
    [ -44.0,   -69.0,   -37.0], \
    [  77.0,   -37.0,   -27.0]]
 
# t4 = t_{ 3/2, sqrt(3)/2 }  (t5 in the paper)
  thop[:,:,3] =  \
   [[ -22.0,    10.0,   -18.0 ], \
    [  10.0,    18.0,    26.0 ], \
    [ -18.0,    26.0,   -11.0 ] ]

# t5 = t_{ 0, sqrt(3) } (t6 in the paper)
  thop[:,:,4] =  \
   [[ -21.0,   -21.0,     0.0 ], \
    [ -21.0,   -26.0,     0.0 ], \
    [   0.0,     0.0,    33.0 ]] 

# t6 = t_{ 3/2, -sqrt(3)/2 } (t4 in the paper)
  thop[:,:,5] =  \
   [[ -22.0,    10.0,    18.0 ], \
    [  10.0,    18.0,   -26.0 ], \
    [  18.0,   -26.0,   -11.0 ] ]

# t7 = t_{ 1, sqrt(3) }
  thop[:,:,6] =  \
   [[ -25.0,   -21.0,   -36.0 ], \
    [ -21.0,   -17.0,   -39.0 ], \
    [ -36.0,   -39.0,   -62.0 ]]

# t8 = t_{ 2, 0 } 
  thop[:,:,7] =  \
   [[ -25.0,    42.0,     0.0 ], \
    [  42.0,   -85.0,     0.0 ], \
    [   0.0,     0.0,     6.0 ]]

# t9 = t_{ 1, -sqrt(3) } = t_{ 2, 0, 0 }
  thop[:,:,8] =  \
   [[ -25.0,   -21.0,    36.0 ], \
    [ -21.0,   -17.0,    39.0 ], \
    [  36.0,    39.0,   -62.0 ]]

# Proper sign for hopping parameters
  thop = -thop / 1000

# *************** DEBUG ************
#  thop[:,:,3:] = 0.0
#  thop[1:,1:,0:3] = 0.0
#  thop[1:,0,0:3] = 0.0
#  thop[0,1:,0:3] = 0.0
#  thop[2:,0:2,0:3] = 0.0
#  thop[0:2,2:,0:3] = 0.0

  eonsite = array( [0.0, -53.0, -53.0] ) / 1000

  return thop, eonsite
 
def Init_Hops_NaxCoO2_Korshunov_L3_x70():
  """
  Initializes hopping parameters for NaxCoO2 with 3 orbitals,
  a1g, eg, eg, included.
  Parameters are from Korshunov et al. (x = 0.7)
  """

  thop = zeros((3,3,9), float)

#
# Notations:
#   1. Plain index notation (like in the notes).
#   2. Absolute normalized coordinates.
#
#      (The data is from Korshunov et al.).
#
# t1 = t_{ 1/2, sqrt(3)/2 }
  thop[:,:,0] = \
   [[ 105.0,   -52.0,   -90.0], \
    [ -52.0,   -68.0,    16.0], \
    [ -90.0,    16.0,   -48.0]]

# t2 = t_{ 1, 0 }
  thop[:,:,1] =  \
   [[ 105.0,   104.0,     0.0 ], \
    [ 104.0,   -39.0,     0.0 ], \
    [   0.0,     0.0,   -77.0 ] ]

# t3 = t_{ 1/2, -sqrt(3)/2 } (the sign of 90.0 seems to be wrong in the paper)
  thop[:,:,2] =  \
   [[ 105.0,   -52.0,    90.0], \
    [ -52.0,   -68.0,   -16.0], \
    [  90.0,   -16.0,   -49.0]]
 
# t4 = t_{ 3/2, sqrt(3)/2 }  (t5 in the paper)
  thop[:,:,3] =  \
   [[ -18.0,     7.0,   -13.0 ], \
    [   7.0,    14.0,    20.0 ], \
    [ -13.0,    20.0,    -9.0 ] ]

# t5 = t_{ 0, sqrt(3) } (t6 in the paper)
  thop[:,:,4] =  \
   [[ -18.0,   -15.0,     0.0 ], \
    [ -15.0,   -20.0,     0.0 ], \
    [   0.0,     0.0,    26.0 ]] 

# t6 = t_{ 3/2, -sqrt(3)/2 } (t4 in the paper)
  thop[:,:,5] =  \
   [[ -18.0,     7.0,    13.0 ], \
    [   7.0,    14.0,   -20.0 ], \
    [  13.0,   -20.0,    -9.0 ] ]

# t7 = t_{ 1, sqrt(3) }
  thop[:,:,6] =  \
   [[ -22.0,   -16.0,   -28.0 ], \
    [ -16.0,   -15.0,   -34.0 ], \
    [ -28.0,   -34.0,   -54.0 ]]

# t8 = t_{ 2, 0 } 
  thop[:,:,7] =  \
   [[ -22.0,    33.0,     0.0 ], \
    [  33.0,   -73.0,     0.0 ], \
    [   0.0,     0.0,     5.0 ]]

# t9 = t_{ 1, -sqrt(3) } = t_{ 2, 0, 0 }
  thop[:,:,8] =  \
   [[ -22.0,   -16.0,    28.0 ], \
    [ -16.0,   -15.0,    34.0 ], \
    [  28.0,    34.0,   -54.0 ]]

# Proper sign for hopping parameters
  thop = -thop / 1000

# *************** DEBUG ************
#  thop[:,:,6:] = 0.0
#  thop[1:,1:,0:3] = 0.0
#  thop[1:,0,0:3] = 0.0
#  thop[0,1:,0:3] = 0.0
#  thop[2:,0:2,0:3] = 0.0
#  thop[0:2,2:,0:3] = 0.0

  eonsite = array( [0.0, -19.0, -19.0] ) / 1000

  return thop, eonsite
 
# ********************************************************************
# ********************************************************************
#
# class Band_Structure
#
# ********************************************************************
# ********************************************************************
class Band_Structure:
  """
  Generic class for band structure construction and manipulation.
  """

# ********************************************************************
#
# __init__
#
# ********************************************************************
  def __init__( self ):

    self.delta = 0.01

    pass


# ********************************************************************
#
# Eigen
#
# ********************************************************************
  def Eigen( self, k ):
    """
    Diagonalizes the Hamiltonian at a given point 'k'
    """

    ham = self.Hamiltonian( k )

# This small correction is needed to avoid abrupt changes
# in the eigenvectors when passing through the k-point
# with H(k) = 0
    if sum( abs( ham[:, :] ) ) < 1e-20:
      ham += 1e-20

    [eigv, eigvecs] = linalg.eigh( ham )
  
    ek = real( eigv )

    return ek, eigvecs


# ********************************************************************
#
# Spec_Fun
#
# ********************************************************************
  def Spec_Fun( self, Chem_pot, k, en, Sigma_k_fun = None, unfold = True, 
                orbs = None, site_block = True ):
    """
    Evaluates the spectral function for a given k-point 'k'
    and energy mesh 'en'.

       * unfold = True : apply the up-folding procedure
         unfold = False: do not apply up-folding

    """

#
# Save the coordinates of the original k-point
#
    k1 = k[0]
    k2 = k[1]

    k1o = k1
    k2o = k2
#  
#   NBZ
    b1 = self.BZ_bas[0, :] / 2 / pi
    b2 = self.BZ_bas[1, :] / 2 / pi
#
#   SBZ
    B = self.SBZ_bas

#   Find K = [k], where k = 1/(2\pi)*(k1*b1 + k2*b2)

    kxy = k1*b1 + k2*b2
    
    mshift = self.Shift_to_BZ(kxy, B)
#  
#   Transform k1,k2 into a point inside the SBZ
#  
    k1 += pi*mshift[0]
    k2 += pi*mshift[1]

    kv = array( [k1, k2] )
#  
#  
    ek, eigvecs = self.Eigen( kv )

    if not Sigma_k_fun is None:
      Sigma = Sigma_k_fun( en, kv )
    else:
      Sigma = None

#  
#   Matrix elements for un-folding
#  
    nm = self.nb * self.no

    if unfold:
      matel = zeros( (nm, self.no), dtype = complex)

      for i in range( self.no ):
        for ij in range( self.nb * self.no ):
          for il in range( self.nb ):
            expf = exp( 1.0j * (k1o * self.sbas[il,0] + k2o * self.sbas[il,1]))
            if site_block:
# Site-blocks order (block corresponds to a site, element within the block -- to an orbital)
              ibl = il * self.no + i
              matel[ij, i] += eigvecs[ibl, ij]*expf
# Orbitals-blocks order (block corresponds to an orbital, element within the block -- to a site)
            else:
              ibl = il + i * self.nb
              matel[ij, i] += eigvecs[ibl, ij]*expf

      matel /= 2.0

    else:
      matel = eye( nm, dtype = complex)

   
    unit = eye( nm, dtype = float )
    ek_mat = diag( ek[:] - Chem_pot, 0 )

    if not Sigma is None:
      delta = 0.0

      SigmaTrans = lambda ie: dot( dot( eigvecs.T.conj(), Sigma[ie,:,:] ), eigvecs )
      Sigma_bloch = array(map(SigmaTrans, arange(len(en))))

      if unfold:
#        G_Fun = vectorize(lambda ie: -dot( matel.T.conj(), dot ( diag( diag( 
#                          linalg.inv( (en[ie] + 1.0j*self.delta )*unit - ek_mat -
#                  Sigma_bloch[ie,:,:] ) ), 0 ), matel ) ), otypes = [ndarray])
        G_Fun = vectorize(lambda ie: -dot( matel.T.conj(), dot ( ( 
                          linalg.inv( (en[ie] + 1.0j*self.delta )*unit - ek_mat -
                  Sigma_bloch[ie,:,:] ) ), matel ) ), otypes = [ndarray])
      else:
        G_Fun = vectorize( lambda ie: 
                         -linalg.inv( (en[ie] + 1.0j*self.delta )*unit - ek_mat -
                  Sigma_bloch[ie,:,:]), otypes = [ndarray] )

    else:
# Sigma is None
      if unfold:
        G_Fun = vectorize(lambda ie: -dot( matel.T.conj(), dot ( 
                          linalg.inv( (en[ie] + 1.0j*self.delta )*unit - ek_mat
                   ), matel ) ), otypes = [ndarray])
#        G_Fun = vectorize(lambda ie: -dot( matel.T.conj(), dot ( diag( diag( 
#                          linalg.inv( (en[ie] + 1.0j*self.delta )*unit - ek_mat
#                   ) ), 0 ), matel ) ), otypes = [ndarray])
      else:
        G_Fun = vectorize( lambda ie: 
                         -linalg.inv( (en[ie] + 1.0j*self.delta )*unit - ek_mat), 
                          otypes = [ndarray] )


    Akw = zeros(shape(en), dtype = float64)

    temp = G_Fun( arange(len(en)) )

    temp2 = zeros( (len(en), ) + temp[0].shape, dtype = complex )

    if orbs is None:
      nto = temp[0].shape[0]
      orbs = arange(nto)
    else:
      orbs = orbs if isinstance(orbs, ndarray) else array(orbs)
      
    for ie in range(len(en)):
      temp2[ie,:] = temp[ie][:]
      Akw[ie] = imag( sum( temp2[ie, orbs, orbs] ) )
#      Akw[ie] = imag( sum( diag( temp2[ie, :, :] ) ) )
#      Akw[ie] = imag( sum( diag( temp2[ie, 1:, 1:] ) ) )
#      Akw[ie] = imag( ( temp2[ie, 0, 0] ) )

    return Akw

# ********************************************************************
#
# Shift_to_BZ
#
# ********************************************************************
  def Shift_to_BZ( self, k, b):
    """
    Shifts an arbitrary k-vector 'k' to the BZ defined
    by 'b'

      * k: general k-vector [kx,ky]
      * b: reciprocal lattice vectors

    Returns:
      * integer numbers m1, m2 such that
        a vector m1*b1 + m2*b2 shifts k into the BZ(b1,b2)

    """

    b1 = b[:,0]
    b2 = b[:,1]

    d = array([dot(k,b1), dot(k,b2)])

    a11 = dot(b1,b1)
    a12 = dot(b1,b2)
    a22 = dot(b2,b2)

    a = array([[a11, a12], [a12, a22]])

    [k1, k2] = linalg.solve(a,d)

    m1 = - int(k1)
    m2 = - int(k2)

    return array( [m1, m2] )

# ********************************************************************
#
# HMLT_Save
#
# ********************************************************************
  def HMLT_Save( self, base, nk ):
    """
    Saves a Hamiltonian in HMLT format
      HMLT.<base> - k-points and upper triangular part of the Hamiltonian
    """

    Fname = "HMLT." + base

    print "  Storing the Hamiltonian to file:", Fname

    fhmlt = open(Fname, "wt")

# Generate k-mesh and calculate Hamiltonian, eigenvalues and
# eigenvectors.
    self.Init_Hamiltonian( nk )

    kweight = 1.0 / self.nktot

    nL = self.nb * self.no
# Initialize files:
#    Nktot, Ntet, Ns, Norb, Efermi
    fhmlt.write("%10d  %8d\n%10d  %8d\n%24.15e\n"%(self.nktot, 0, 1, nL, 0.0))

#
# Store the data to the files.
#
# Weights
    for ik in range( self.nktot ):
      fhmlt.write("%24.15e"%( kweight ) )
      if ((ik + 1) / 3) * 3 == ik + 1:
        fhmlt.write("\n")

    if (self.nktot / 3) * 3 != self.nktot:
        fhmlt.write("\n")

#
# K-points
#
    for ik2 in range(nk):
      for ik1 in range(nk):
        kv = self.kmesh[ik1, ik2, :]

        fhmlt.write("%24.15e%24.15e%24.15e\n"%(kv[0], kv[1], kv[2]))

    iv = 0
    for ik2 in range(nk):
      for ik1 in range(nk):
        U = self.Uk[ik1, ik2, :, :]

        for i in range(nL):
          for j in range(i, nL):
            ham = self.hamk[ik1, ik2, i, j]
            fhmlt.write("%24.15e"%( real(ham) ) )
            iv += 1
            if (iv / 3) * 3 == iv:
              fhmlt.write("\n")

        for i in range(nL):
          for j in range(i, nL):
            ham = self.hamk[ik1, ik2, i, j]
            fhmlt.write("%24.15e"%( imag(ham) ) )
            iv += 1
            if (iv / 3) * 3 == iv:
              fhmlt.write("\n")

# ********************************************************************
#
# Generate_Save
#
# ********************************************************************
  def Generate_Save( self, base, nk ):
    """
    Generate and store data into files:
      <base>.eig - eigenvalues
      <base>.plo - projectors
      <base>.kpoints - k-points and weights
      <base>.matel - matrix elements
    """

    def transf_apply( A, U ):
      return dot( U.conj().T, dot( A, U ) )

    FeigName = base + ".eig"
    FmatelName = base + ".matel"
    FploName = base + ".plo"
    FkptsName = base + ".kpoints"

    feig = open(FeigName, "wt")
    fmatel = open(FmatelName, "wt")
    fplo = open(FploName, "wt")
    fkpts = open(FkptsName, "wt")

# Generate k-mesh and calculate Hamiltonian, eigenvalues and
# eigenvectors.
    self.Init_Hamiltonian( nk )

    kweight = 1.0 / self.nktot
    nL = self.nb * self.no

# Initialize files
    fplo.write("%10d  %8d  %8d  %8d  %8d\n\n\n"%(nL, nL, self.nktot, 1, 1))

# Store the data to the files.
# Projectors are calculated as eigenvector matrices.
    ikp = 0
    for ik2 in range(nk):
      for ik1 in range(nk):
        ikp += 1

        kv = self.kmesh[ik1, ik2, :]

        fkpts.write("%20.15f %20.15f %20.15f %20.15f\n"%(kv[0], kv[1], kv[2], kweight))

        for ib in range(self.nb * self.no):
          feig.write("%10d %9d %29.20f\n"%( ib+1, ikp, self.eigk[ik1, ik2, ib]) )

        fplo.write("\n")

        for io in range( self.Uk.shape[2] ):
          for ib in range( self.Uk.shape[3] ):
            proj = self.Uk[ik1, ik2, io, ib]
            fplo.write("%30.20f %30.20f\n"%(real(proj), imag(proj)))
#
# Matrix elements
#
    matel = zeros( (2, nL, nL), dtype = complex )

    dk = 2 * pi / ( nk  )

    melk = self.Numerical_Matrix_Elements( [dk, dk] )

    print
    print "  Storing matrix elements..."

    ikp = 0
    averr = 0.0
    for ik2 in range(nk):
      for ik1 in range(nk):
#
# Debug code left here for testing matrix elements for new geometries
#        kv = self.kmesh[ik1, ik2, :]
#        [vx, vy] = [ 0.4 * sin( kv[0] ), 0.4 * sin( kv[1] ) ]
#
#        err = (vx - melk[0, ik1, ik2, 0, 0])**2 + (vy - melk[1, ik1, ik2, 0, 0])**2
#        averr += err
#        assert err < 1.0e-3/nk, \
#               "  Something is wrong: %d  %d"%(ik1, ik2)

        ikp += 1

        for i in range(2):
          matel[i, :, :] = transf_apply( melk[i, ik1, ik2, :, :], self.Uk[ik1, ik2, :, :] )

        for ib2 in range(nL):
          for ib1 in range(nL):
            for i in range(3):
              if i < 2:
                vr = real(matel[i, ib1, ib2])
                vi = imag(matel[i, ib1, ib2])
              else:
                vr = 0.0
                vi = 0.0

              fmatel.write("%10d %9d %9d %9d %29.20f %29.20f\n"\
                           %(ib1+1, ib2+1, ikp, i+1, vr, vi))

        fmatel.write("\n")

    feig.close()
    fmatel.close()
    fkpts.close()
    fplo.close()


# ********************************************************************
#
# Numerical_Matrix_Elements
#
# ********************************************************************
  def Numerical_Matrix_Elements( self, dk ):
    """
    Calculate matrix elements by taking a numerical derivative
    of the Hamiltonian plus cluster term.
    """

# Coefficients for the 3d order formula for the numerical derivative
    coefs = array( [1, -8, 8, -1], dtype = float) / 12.0

    abas = linalg.inv( self.SBZ_bas )

    nkx = self.hamk.shape[0]
    nky = self.hamk.shape[1]
    nb = self.hamk.shape[2]

    A = array( [abas[:,0] * nkx, abas[0:,1] * nky] ).T

    velk = zeros( (2, nkx, nky, nb, nb), dtype = complex )
    sij = zeros( (2, nb, nb), dtype = float)

    for ib2 in range(nb):
      for ib1 in range(nb):
        sij[:, ib1, ib2] = self.svec[:, ib1] - self.svec[:, ib2]

    p1 = zeros( 4, dtype = int )
    p2 = zeros( 4, dtype = int )

    kscal_1 = 2 * pi / nkx / dk[0]
    kscal_2 = 2 * pi / nky / dk[1]

    for ik2 in range(nky):
      for ik1 in range(nkx):
        p1[0] = (ik1 - 2 + nkx) % nkx
        p1[1] = (ik1 - 1 + nkx) % nkx
        p1[2] = (ik1 + 1) % nkx
        p1[3] = (ik1 + 2) % nkx

        p2[0] = (ik2 - 2 + nky) % nky
        p2[1] = (ik2 - 1 + nky) % nky
        p2[2] = (ik2 + 1) % nky
        p2[3] = (ik2 + 2) % nky

        for ib2 in range(nb):
          for ib1 in range(nb):
            v1 = sum( self.hamk[p1, ik2, ib1, ib2] * coefs ) * kscal_1

            v2 = sum( self.hamk[ik1, p2, ib1, ib2] * coefs ) * kscal_2

            velk[0, ik1, ik2, ib1, ib2] = A[0,0] * v1 + A[0,1] * v2
            velk[1, ik1, ik2, ib1, ib2] = A[1,0] * v1 + A[1,1] * v2

        velk[0, ik1, ik2, :, :] += 1.0j * sij[0, :, :] * self.hamk[ik1, ik2, :, :]
        velk[1, ik1, ik2, :, :] += 1.0j * sij[1, :, :] * self.hamk[ik1, ik2, :, :]
    
    return velk

# ********************************************************************
#
# Read_External
#
# ********************************************************************
  def Read_External( self, base, nk, BZ_bas ):
    """
    Reads externally generated band-structure data.
    Uniform 2D k-mesh with 'nk' points in each direction is assumed.

    The data is read from files:
      <base>.eig - eigenvalues
      <base>.plo - projectors
      <base>.kpoints - k-points and weights
      <base>.matel - matrix elements
    """

    def transf_apply( A, U ):
      return dot( U.conj().T, dot( A, U ) )

# Convert a string to a float and to an int
    tofloat = lambda s: float(s)
    toint = lambda s: int(s)

# Convert a list of two strings to a complex
    tocomplex = lambda s: complex(float(s[0]), float(s[1]))

    FeigName = base + ".eig"
    FmatelName = base + ".matel"
    FploName = base + ".plo"
    FkptsName = base + ".kpoints"

    feig = open(FeigName, "rt")
    fmatel = open(FmatelName, "rt")
    fplo = open(FploName, "rt")
    fkpts = open(FkptsName, "rt")

    line = fplo.readline()

# Read dimensions from the first line of .plo-file
    [nb, norb, nktot, nspin, ndim] = map(toint, line.split())

    print
    print "Dims:", nb, norb, nktot, nspin, ndim

    wks = zeros(nktot, dtype = float)
    kpts = zeros((3, nktot), dtype = float)
    plo = zeros((nktot, norb, nb), dtype = complex)
    eigs = zeros((nktot, nb), dtype = float)
    fvel = zeros((3, nktot, nb, nb), dtype = complex)

# Read two empty lines
    line = fplo.readline()
    line = fplo.readline()

    for ik in range(nktot):
      for ib in range(nb):
        sline = feig.readline().split()
        eigs[ik, ib] = float(sline[2])

# Read trailing empty line
      line = fplo.readline()

      for iorb in range(norb):
        for ib in range(nb):
          sline = fplo.readline().split()
          plo[ik, iorb, ib] = tocomplex(sline)

      sline = fkpts.readline().split()
      wks[ik] = float(sline[3])
#    kpts_1d[ik] = float(ik) / (nktot - 1)
      kpts[0, ik] = float(sline[0])
      kpts[1, ik] = float(sline[1])
      kpts[2, ik] = float(sline[2])

#
# Read matrix elements
#
    for ik in range(nktot):
      for ib2 in range(nb):
        for ib1 in range(nb):
          for i in range(3):
            sline = fmatel.readline().split()
            fvel[i, ik, ib1, ib2] = tocomplex(sline[4:6])

      fmatel.readline()

    feig.close()
    fmatel.close()
    fkpts.close()
    fplo.close()

    print
    print "  Checking orthogonality of projectors..."

    averr = 0.0
    unit = eye( nb, dtype = complex )
    for ik in range(nktot):
      plok = plo[ik, :, :]
      prod = dot( plok.conj().T, plok )
      err = sqrt( sum( abs(prod - unit)**2 ) )
      averr += err
      
    averr /= nktot

    print "  Average orthogonality error:", averr

#
# Recover the topology of the k-mesh
#
    self.kmesh = zeros( (nk, nk, 3), dtype = float )
    self.hamk = zeros( (nk, nk, norb, norb), dtype = complex )
    self.eigk = zeros( (nk, nk, nb), dtype = float )

    dk = 1.0 / ( nk  )
    kr1 = linspace(-0.5, 0.5 - dk, nk)
    kr2 = linspace(-0.5, 0.5 - dk, nk)

    for ik2 in range(nk):
      for ik1 in range(nk):

        k1 = kr1[ik1]
        k2 = kr2[ik2]

        self.kmesh[ik1, ik2, 0:2] = BZ_bas[0, 0:2] * k1 + BZ_bas[1, 0:2] * k2

# Create mapping between the k-point list and k-mesh
    print
    print "  Constructing the KD-tree..."
            
    kdtree = KDTree( kpts[0:2,:].T )

    self.kmap = -1 * ones( (nk, nk), dtype = int )

    for ik2 in range(nk):
      for ik1 in range(nk):
        k = self.kmesh[ik1, ik2, 0:2]
        d, i = kdtree.query(k)
        if d < 1e-8:
          self.kmap[ik1, ik2] = i
        else:
          for nx in range(-1, 2):
            for ny in range(-1, 2):
              ks = k + nx * BZ_bas[0, 0:2] + ny * BZ_bas[1, 0:2]
              d, i = kdtree.query(ks)
              if d < 1e-8:
                self.kmap[ik1, ik2] = i

        assert self.kmap[ik1, ik2] >= 0, "  Cannot find point (%f, %f)"%(k[0], k[1])

#    print "  Mapping:"
#    for ik2 in range(nk):
#      for ik1 in range(nk):
#        i = self.kmap[ik1, ik2]
#        print ik1, ik2, self.kmesh[ik1, ik2, 0:2], "  --> ", i, kpts[0:2, i]

# ********************************************************************
#
# Init_Hamiltonian
#
# ********************************************************************
  def Init_Hamiltonian( self, nk ):
    """
    Generate the k-mesh and Hamiltonian.
      * nk -- number of k-points in each direction
    """

    self.nk = nk

    nb = self.nb
    no = self.no
    nL = nb * no

    dk = 1.0 / ( nk  )
    kr1 = linspace(-0.5, 0.5 - dk, nk)
    kr2 = linspace(-0.5, 0.5 - dk, nk)

    b1 = self.SBZ_bas[0, :] / 2 / pi
    b2 = self.SBZ_bas[1, :] / 2 / pi

    print "b1,b2:", b1,b2
    self.nktot = nk*nk

    self.hamk = zeros( (nk, nk, nL, nL), dtype = complex )
    self.eigk = zeros( (nk, nk, nL), dtype = float )
    self.kmesh = zeros( (nk, nk, 3), dtype = float )
    self.Uk = zeros( (nk, nk, nL, nL), dtype = complex )


    for ik2 in range(nk):
      for ik1 in range(nk):

        k1 = 2 * pi * kr1[ik1]
        k2 = 2 * pi * kr2[ik2]
        kp = array( [k1, k2] )

        self.kmesh[ik1, ik2, 0:2] = b1[0:2] * k1 + b2[0:2] * k2
        self.hamk[ik1, ik2, :, :] = self.Hamiltonian( kp )

        eigv, eigvecs = self.Eigen( kp )
        self.eigk[ik1, ik2, :] = eigv
        self.Uk[ik1, ik2, :, :] = eigvecs

##################################################################
#
# Crystal_Field
#
##################################################################
  def Crystal_Field(self):
    """
    Crystal field calculation.
    """

    nk = self.nk

    kweight = 1.0 / self.nktot

    ndim = bs_kag.no * bs_kag.nb

    hdiag = zeros( ndim, complex )

    for ik1 in range(nk):
      for ik2 in range(nk):
        ham = self.hamk[ik1, ik2, :, :]
        hdiag += diag( ham )

    hdiag *= kweight

    print "  Diagonal part of H(R = 0):"
    print hdiag

    return hdiag

##################################################################
#
# Crystal_Field
#
##################################################################
  def Density_of_States(self, Efermi):
    """
    The density of states.
    """

    assert 0, "  The method is not yet implemented"

#
# Density and number of states
#
#  bs_tri = BS_Tri4_L3( V4 )

    bs_kag = BS_Kag( )

    nb = 3

    V = 0.3
    Sigma_fun = Read_Sigma("kag0.67/V0.00/SigmaRe_site", nb)

    print " Im Sigma( 0.0 ):"
    print imag( Sigma_fun(0.0) )

    def Sigma_k_fun( en, kv ):
#
# Periodized Sigma
#
      k1 = kv[0]
      k2 = kv[1]
      k21 = k2 - k1

      Sigma = Sigma_fun(en)
      Sigma_k = array( Sigma )

# Self-energy periodization
#    Sig_0 = Sigma[:, 0, 0] - Efermi
#    Sig_1 = Sigma[:, 0, 1]
#
#    Sigma_k[:, 0, 0] = Sig_0
#    Sigma_k[:, 1, 1] = Sig_0
#    Sigma_k[:, 2, 2] = Sig_0
#    Sigma_k[:, 1, 0] = 0.5 * Sig_1 * (1.0 + exp( -2.0j * k1 ) )
#    Sigma_k[:, 2, 0] = 0.5 * Sig_1 * (1.0 + exp( -2.0j * k2 ) )
#    Sigma_k[:, 2, 1] = 0.5 * Sig_1 * (1.0 + exp( -2.0j * k21 ) )
#    Sigma_k[:, 0, 1] = 0.5 * Sig_1 * (1.0 + exp(  2.0j * k1 ) )
#    Sigma_k[:, 0, 2] = 0.5 * Sig_1 * (1.0 + exp(  2.0j * k2 ) )
#    Sigma_k[:, 1, 2] = 0.5 * Sig_1 * (1.0 + exp(  2.0j * k21 ) )


# Cumulant periodization
      M = zeros( Sigma.shape, dtype = complex )
      M_k = zeros( Sigma.shape, dtype = complex )
      umat = eye( nb, dtype = complex )

      for ie, e in enumerate( en ):
        M[ie, : :] = linalg.inv( (e + Efermi) * umat  - Sigma[ie, :, :] )

      M_0 = M[:, 0, 0]
      M_1 = M[:, 0, 1]

      M_k[:, 0, 0] = M_0
      M_k[:, 1, 1] = M_0
      M_k[:, 2, 2] = M_0
      M_k[:, 1, 0] = 0.5 * M_1 * (1.0 + exp( -2.0j * k1 ) )
      M_k[:, 2, 0] = 0.5 * M_1 * (1.0 + exp( -2.0j * k2 ) )
      M_k[:, 2, 1] = 0.5 * M_1 * (1.0 + exp( -2.0j * k21 ) )
      M_k[:, 0, 1] = 0.5 * M_1 * (1.0 + exp(  2.0j * k1 ) )
      M_k[:, 0, 2] = 0.5 * M_1 * (1.0 + exp(  2.0j * k2 ) )
      M_k[:, 1, 2] = 0.5 * M_1 * (1.0 + exp(  2.0j * k21 ) )

      for ie, e in enumerate( en ):
# N.B. Efermi is subtracted!
        Sigma_k[ie, : :] = e - linalg.inv( M_k[ie, :, :] )

      return Sigma_k

    nk = 101
    nkz = 1
    ne = 501

    emin = -1.5
    emax =  1.5

    kr1 = linspace(-pi, pi, nk)
    kr2 = linspace(-pi, pi, nk)

    ens = linspace(emin, emax, ne)  
    deltE = ens[2] - ens[1]
#
#   
#  Sigma = zeros((len(ens), nb), dtype = complex)

    nktot = nk*nk*nkz

    kweight = 1.0/nktot
#
# 2 for spin DoF
# 1 / 4 for 4-site cluster
#
    fac = 2.0*kweight/pi/nb

    dos = zeros(ne, dtype = float64)
    nos = zeros(ne, dtype = float64)

    if nkz > 1:
      nkz2 = nkz / 2
    else:
      nkz2 = 1

    for ik1 in range(nk):
      for ik2 in range(nk):
        
        k1 = kr1[ik1]
        k2 = kr2[ik2]

        for ikz in range(nkz):
          kz = float(ikz - nkz/2) / ( nkz2 ) * pi
          ekz = -2 * thopz * cos(kz) * 0.0

          ensz = ens + ekz
          
#        Akk = bs_tri.Spec_Fun( array( [k1, k2] ), ensz, upfold = 0 )
          Akk = bs_kag.Spec_Fun( array( [k1, k2] ), ensz, upfold = 0 )
#        Akk = bs_kag.Spec_Fun( array( [k1, k2] ), ensz, Sigma_k_fun, upfold = 0 )
#        Akk = tri4_sf(k1,k2,ensz,None)
#        Akk = tri_sf(k1,k2,ensz,None)
          dos += fac*Akk

    for ie in range(1,ne):
      nos[ie] = nos[ie-1] + (dos[ie-1] + dos[ie])/2

    nos *= deltE

    print "The total number of states: ", sum(dos)*deltE, nos[ne-1]

#  fo = open("dos_kag0.56_V%.2f.dos"%(V), "wt")
    fo = open("dos_kag_t3_bare.dos", "wt")

    fo.write("# DOS and NOS for Kagome lattice with V = %.2f\n"%(V))

    for ie in range(ne):
      fo.write("%.5f   %.10f   %.10f\n"%(ens[ie], dos[ie], nos[ie]))

    fo.close()
    d = Gnuplot.Data(ens, dos, with_ = "l")
    d2 = Gnuplot.Data(ens, nos, with_ = "l")

    graph.plot(d, d2)
#
#


##################################################################
#
# set_kvec
#
##################################################################
  def set_kvec(self, ki, kf, nk):
    """
    Returns an array of 2D vectors with nk vectors ranging from ki to kf.
    """

    b1 = self.BZ_bas[0, :] / 2 / pi
    b2 = self.BZ_bas[1, :] / 2 / pi

    kxi = ki[0]*b1[0] + ki[1]*b2[0]
    kyi = ki[0]*b1[1] + ki[1]*b2[1]

    dk1 = (kf[0] - ki[0])/(nk - 1)
    dk2 = (kf[1] - ki[1])/(nk - 1)

    kvec = zeros((nk,2), float64)
    kmod = zeros(nk, float64)

    for i in range(nk):
      k1 = ki[0] + dk1*i
      k2 = ki[1] + dk2*i

      kx = k1*b1[0] + k2*b2[0]
      ky = k1*b1[1] + k2*b2[1]

      kdist = sqrt((kx - kxi)**2 + (ky - kyi)**2)

      kmod[i] = kdist
      kvec[i][0] = k1
      kvec[i][1] = k2

    return kmod, kvec


##################################################################
#
# SF_kpoint
#
##################################################################
  def SF_kpoint(self, Chem_pot, kp, ens, Sigma_k_fun, plotter, 
                unfold=True, orbs=None):
    """
    Plot the spectral function at a particular k-point 'kp'.
    """
#
# Calculate eigenvalues at the corresponding k-point
#
    ne = ens.shape[0]

    Asf = zeros(ne)
#
# Spectral function in k-e planes
#
    Asf[:] = self.Spec_Fun(Chem_pot, array( [kp[0], kp[1]] ), ens, Sigma_k_fun, 
                           unfold, orbs)

# Filter "bad" points (usually appearing due to non-causality of Sigma(k) )
    Asf[ Asf < 0.0 ] = 0.0

    plotter.Plot_SF_kpoint(ens, Asf, 0.0)

##################################################################
#
# SF_band_structure
#
##################################################################
  def SF_band_structure(self, Chem_pot, nk, ens, Sigma_k_fun, plotter,
                        unfold=True, orbs=None, Stacks3D=False):
    """
    Plot the band structure in k-e plane.
    """

    ne = ens.shape[0]

    klen, kdat, nr, lett = self.High_Symmetry_kPoints(nk)

    Asf = zeros((ne, nk, nr), float64)

#
# Spectral function in the k-e plane
#
    ip = 0

    for ir in range(nr):
      for ik in range(nk):
        k1 = kdat[ip][0]
        k2 = kdat[ip][1]
        ip += 1
        
        Asf[:,ik,ir] = self.Spec_Fun( Chem_pot, array( [k1, k2] ), ens, 
                         Sigma_k_fun, unfold=unfold, orbs=orbs )

    Asf[ Asf < 0.0 ] = 0.0

    if Stacks3D:
      plotter.Plot_SF_stack(klen, ens, Asf, 0.0, nk, nr, lett)
    else:
      plotter.Plot_SF(klen, ens, Asf, 0.0, nk, nr, lett)

##################################################################
#
# SF_Fermi_surface
#
##################################################################
  def SF_Fermi_surface(self, Chem_pot, nk, kmax, en_fermi, Sigma_k_fun, 
                       plotter, unfold=True, orbs=None):
    """
    Plot the spectral function in the k-k plane (Fermi surface).
      * en_fermi -- a small set of energy points around the Fermi level
                    over which the output is averaged; needed for smoothing.
    """

    krx = linspace(-kmax, kmax, nk)
    kry = linspace(-kmax, kmax, nk)

    Afermi = zeros((nk, nk), dtype = float64)

    abas = 2 * pi * linalg.inv(self.BZ_bas)
    print "  abas:", abas

    for ik1 in range(nk):
      for ik2 in range(nk):
#
# Find k1, k2 for a point (krx[ik1], kry[ik2])
#
        d = array([krx[ik1], kry[ik2]])
        [k1, k2] = dot(abas.T, d)

        Akk = self.Spec_Fun(Chem_pot, array([k1, k2]), en_fermi, Sigma_k_fun, 
                            unfold=unfold, orbs=orbs)
        Afermi[ik1,ik2] = mean(Akk)
      
    plotter.Plot_Fermi_surface(krx, kry, Afermi)


# ********************************************************************
# ********************************************************************
#
# class BS_Tri4
#
# ********************************************************************
# ********************************************************************
class BS_Tri4( Band_Structure ):
  """
  4-site rhomb cluster on the triangular lattice
  """

# ********************************************************************
#
# __init__
#
# ********************************************************************
  def __init__( self, V4 ):
    """
    """
    self.nb = 4
    self.no = 1

    self.thop = -0.134
    self.thop2 = 0.028
    self.thop3 = 0.024
    self.thopz = 0.015

    self.V4 = V4
#
# BZ basis vectors
#
    self.BZ_bas = 2 * pi * array( [[1, -1/sqrt(3)], \
                                   [0,  2/sqrt(3)]] )

#
# Super BZ (smaller BZ of the supercell)
#
    self.SBZ_bas = pi * array( [[1, -1/sqrt(3)], \
                                [0,  2/sqrt(3)]] )
# 
#  Cluster basis vectors 
#  (defined in terms of a1, a2 - translation vectors of the original lattice)
# 
    self.sbas = array([[0.0, 0.0],
                       [1.0, 0.0],
                       [0.0, 1.0],
                       [1.0, 1.0]])

    self.svec = zeros( (2, 4), dtype=float )

    Abas = 2 * pi * linalg.inv( self.BZ_bas )

    print "  Translation vectors of the underlying lattice:"
    print Abas

    for i in range(4):
      self.svec[:, i] = dot( Abas, self.sbas[i, :] )

    Band_Structure.__init__( self )

# ********************************************************************
#
# Hamiltonian
#
# ********************************************************************
  def Hamiltonian( self, kp ):
    """
    Evaluates the Hamiltonian for a given point 'k'
    """

    k1 = kp[0]
    k2 = kp[1]

    emat = zeros((4,4), complex64)
    ematp = zeros((4,4), complex64)
    emat3 = zeros((4,4), complex64)

    emat[0,1] = -( 1 + exp( 2.0j*k1) )
    emat[1,0] = -( 1 + exp(-2.0j*k1) )
    emat[0,2] = -( 1 + exp( 2.0j*k2) )
    emat[2,0] = -( 1 + exp(-2.0j*k2) )
    emat[0,3] = -( exp( 2.0j*k1) + exp( 2.0j*k2) )
    emat[3,0] = -( exp(-2.0j*k1) + exp(-2.0j*k2) )
    emat[1,2] = -( 1 + exp( 2.0j*(k2-k1)) )
    emat[2,1] = -( 1 + exp(-2.0j*(k2-k1)) )
    emat[1,3] = -( 1 + exp( 2.0j*k2) )
    emat[3,1] = -( 1 + exp(-2.0j*k2) )
    emat[2,3] = -( 1 + exp( 2.0j*k1) )
    emat[3,2] = -( 1 + exp(-2.0j*k1) )

    emat *= self.thop

    ematp[1,0] = -( exp(-2.0j*k2) + exp( 2.0j*(k2-k1)) )
    ematp[2,0] = -( exp(-2.0j*k1) + exp(-2.0j*(k2-k1)) )
    ematp[3,0] = -( 1 + exp(-2.0j*(k1+k2)) )
    ematp[2,1] = -( exp(-2.0j*k2) + exp( 2.0j*k1) )
    ematp[3,1] = -( exp(-2.0j*k1) + exp(-2.0j*(k2-k1)) )
    ematp[3,2] = -( exp(-2.0j*k2) + exp( 2.0j*(k2-k1)) )
    ematp[0,1] = conjugate(ematp[1,0])
    ematp[0,2] = conjugate(ematp[2,0])
    ematp[0,3] = conjugate(ematp[3,0])
    ematp[1,2] = conjugate(ematp[2,1])
    ematp[1,3] = conjugate(ematp[3,1])
    ematp[2,3] = conjugate(ematp[3,2])

    ematp *= self.thop2

    e3 = - 2 * self.thop3 * ( cos( 2*k1 ) + cos( 2*k2 ) + cos( 2*(k2 - k1) ) )

    emat3[0,0] = e3
    emat3[1,1] = e3
    emat3[2,2] = e3
    emat3[3,3] = e3

    emat[3,3] = self.V4

    emat_full = emat + ematp + emat3
 
    return emat_full

# ********************************************************************
#
# Spec_Fun
#
# ********************************************************************
  def Spec_Fun( self, Chem_pot, k, en, Sigma_k_fun = None, unfold = True, 
                orbs = None ):
    """
    Evaluates the spectral function for a given k-point 'k'
    and energy mesh 'en'.

       * unfold = True : apply the up-folding procedure
         unfold = False: do not apply up-folding

    """

    return Band_Structure.Spec_Fun(self, Chem_pot, k, en, Sigma_k_fun, 
                                   unfold, orbs, site_block=True)



# ********************************************************************
#
# High_Symmetry_kPoints
#
# ********************************************************************
  def High_Symmetry_kPoints(self, nk):
    """
    Generates lines along high-symmetry directions in k-space.
    """

    nr = 2

    klen = zeros(nk * nr, float64)
    kdat = zeros((nk * nr, 2), float64)

#
# Define 2 k-ranges
#
    ik = 0
    ir = 0
#
# Symbols for special k-points
#
#  lett = ['{/Symbol G}', 'K', 'M', '{/Symbol G}']
    lett = ['M', '{/Symbol G}', 'K']

# M-G
    ki1 = pi
    ki2 = pi
    kf1 = 0.0
    kf2 = 0.0

    i1 = ik
    i2 = ik + nk

    klen[i1:i2], kdat[i1:i2][:] = self.set_kvec([ki1, ki2], [kf1, kf2], nk)

    ir += 1
    ik += nk

# G-K
    ki1 = 0.0
    ki2 = 0.0
    kf1 = 4*pi/3
    kf2 = 2*pi/3
#  kf1 = pi
#  kf2 = pi/2

    i1 = ik
    i2 = ik + nk

    klen[i1:i2], kdat[i1:i2][:] = self.set_kvec([ki1, ki2], [kf1, kf2], nk)

    return klen, kdat, nr, lett


# ********************************************************************
# ********************************************************************
#
# class BS_Tri4_L3
#
# ********************************************************************
# ********************************************************************
class BS_Tri4_L3( Band_Structure ):
  """
  4-site rhomb cluster on the triangular lattice.
  3-oribtals: a1g, eg, eg.
  """

# ********************************************************************
#
# __init__
#
# ********************************************************************
  def __init__( self, V4 ):
    """
    """
    self.nb = 4
    self.no = 3

    self.V4 = V4

#    self.thop, self.eonsite = Init_Hops_NaxCoO2_Korshunov_L3_x70()
    self.thop, self.eonsite = Init_Hops_NaxCoO2_Frank_L3()

#
# BZ basis vectors
#
    self.BZ_bas = 2 * pi * array( [[1, -1/sqrt(3)], \
                                   [0,  2/sqrt(3)]] )

#
# Super BZ (smaller BZ of the supercell)
#
    self.SBZ_bas = pi * array( [[1, -1/sqrt(3)], \
                                [0,  2/sqrt(3)]] )
# 
#  Cluster basis vectors 
#  (defined in terms of a1, a2 - translation vectors of the original lattice)
# 
    self.sbas = array([[0.0, 0.0],
                       [1.0, 0.0],
                       [0.0, 1.0],
                       [1.0, 1.0]])

    Band_Structure.__init__( self )

# ********************************************************************
#
# Hamiltonian
#
# ********************************************************************
  def Hamiltonian( self, kp ):
    """
    Evaluates the Hamiltonian for a given point 'k'
    """

    k1 = kp[0]
    k2 = kp[1]

    thop = self.thop
    eonsite = self.eonsite

    emat0 = zeros((4, 4, 3, 3), complex64)

    emat0[0,0,:] =  diag( eonsite, 0 )
    emat0[1,1,:] =  emat0[0,0,:]
    emat0[2,2,:] =  emat0[0,0,:]
    emat0[3,3,:] =  emat0[0,0,:]

    emat1 = zeros((4, 4, 3, 3), complex64)
    emat2 = zeros((4, 4, 3, 3), complex64)
    emat3 = zeros((4, 4, 3, 3), complex64)

    emat1[1,0,:] = -( thop[:,:,1].T + thop[:,:,1] * exp( -2.0j * k1 ) )
    emat1[2,0,:] = -( thop[:,:,0].T + thop[:,:,0] * exp( -2.0j * k2 ) )
    emat1[3,0,:] = -( thop[:,:,2] * exp( -2.0j * k1 ) + thop[:,:,2].T * exp( -2.0j * k2 ) )
    emat1[2,1,:] = -( thop[:,:,2] + thop[:,:,2].T * exp( -2.0j * (k2-k1) ) )
    emat1[3,1,:] = -( thop[:,:,0].T + thop[:,:,0] * exp( -2.0j * k2 ) )
    emat1[3,2,:] = -( thop[:,:,1].T + thop[:,:,1] * exp( -2.0j * k1 ) )

#    emat1 += emat1.conj().T

    emat2[1,0,:] = -( thop[:,:,4] * exp( -2.0j * k2 ) + thop[:,:,4].T * exp( -2.0j * (k1-k2) ) )
    emat2[2,0,:] = -( thop[:,:,5] * exp( -2.0j * k1 ) + thop[:,:,5].T * exp( -2.0j * (k2-k1) ) )
    emat2[3,0,:] = -( thop[:,:,3].T + thop[:,:,3] * exp( -2.0j * (k1+k2) ) )
    emat2[2,1,:] = -( thop[:,:,3] * exp( -2.0j * k2 ) + thop[:,:,3].T * exp( 2.0j * k1 ) )
    emat2[3,1,:] = -( thop[:,:,5] * exp( -2.0j * k1 ) + thop[:,:,5].T * exp( -2.0j * (k2-k1) ) )
    emat2[3,2,:] = -( thop[:,:,4] * exp( -2.0j * k2 ) + thop[:,:,4].T * exp( -2.0j * (k1-k2) ) )

#    emat2 += emat2.conj().T

    e3 = - ( thop[:,:,7] * exp( -2.0j * k1 ) + \
             thop[:,:,6] * exp( -2.0j * k2 ) + \
             thop[:,:,8] * exp( -2.0j * (k1 - k2) ) )

    emat3[0,0,:] = e3 + e3.conj().T
    emat3[1,1,:] = emat3[0,0,:]
    emat3[2,2,:] = emat3[0,0,:]
    emat3[3,3,:] = emat3[0,0,:]

    emat1[3,3,:] = self.V4 * eye( self.no, dtype = complex64 )
#    emat1[3,3,0,0] = self.V4

    emat0b = bmat( [[emat0[0,0,:], emat0[0,1,:], emat0[0,2,:], emat0[0,3,:]], \
                    [emat0[1,0,:], emat0[1,1,:], emat0[1,2,:], emat0[1,3,:]], \
                    [emat0[2,0,:], emat0[2,1,:], emat0[2,2,:], emat0[2,3,:]], \
                    [emat0[3,0,:], emat0[3,1,:], emat0[3,2,:], emat0[3,3,:]] ] )

    emat1b = bmat( [[emat1[0,0,:], emat1[0,1,:], emat1[0,2,:], emat1[0,3,:]], \
                    [emat1[1,0,:], emat1[1,1,:], emat1[1,2,:], emat1[1,3,:]], \
                    [emat1[2,0,:], emat1[2,1,:], emat1[2,2,:], emat1[2,3,:]], \
                    [emat1[3,0,:], emat1[3,1,:], emat1[3,2,:], emat1[3,3,:]] ] )

    emat2b = bmat( [[emat2[0,0,:], emat2[0,1,:], emat2[0,2,:], emat2[0,3,:]], \
                    [emat2[1,0,:], emat2[1,1,:], emat2[1,2,:], emat2[1,3,:]], \
                    [emat2[2,0,:], emat2[2,1,:], emat2[2,2,:], emat2[2,3,:]], \
                    [emat2[3,0,:], emat2[3,1,:], emat2[3,2,:], emat2[3,3,:]] ] )

    emat3b = bmat( [[emat3[0,0,:], emat3[0,1,:], emat3[0,2,:], emat3[0,3,:]], \
                    [emat3[1,0,:], emat3[1,1,:], emat3[1,2,:], emat3[1,3,:]], \
                    [emat3[2,0,:], emat3[2,1,:], emat3[2,2,:], emat3[2,3,:]], \
                    [emat3[3,0,:], emat3[3,1,:], emat3[3,2,:], emat3[3,3,:]] ] )

    emat1b += emat1b.conj().T
    emat2b += emat2b.conj().T

    emat_full = emat0b + emat1b + emat2b + emat3b
 
    return emat_full

# ********************************************************************
#
# Spec_Fun
#
# ********************************************************************
  def Spec_Fun( self, Chem_pot, k, en, Sigma_k_fun = None, unfold = True, 
                orbs = None ):
    """
    Evaluates the spectral function for a given k-point 'k'
    and energy mesh 'en'.

       * unfold = True : apply the up-folding procedure
         unfold = False: do not apply up-folding

    """

    return Band_Structure.Spec_Fun(self, Chem_pot, k, en, Sigma_k_fun, 
                                   unfold, orbs, site_block=True)


# ********************************************************************
#
# High_Symmetry_kPoints
#
# ********************************************************************
  def High_Symmetry_kPoints(self, nk):
    """
    Generates lines along high-symmetry directions in k-space.
    """

    nr = 2

    klen = zeros(nk * nr, float64)
    kdat = zeros((nk * nr, 2), float64)

#
# Define 2 k-ranges
#
    ik = 0
    ir = 0
#
# Symbols for special k-points
#
#  lett = ['{/Symbol G}', 'K', 'M', '{/Symbol G}']
    lett = ['M', '{/Symbol G}', 'K']

# M-G
    ki1 = pi
    ki2 = pi
    kf1 = 0.0
    kf2 = 0.0

    i1 = ik
    i2 = ik + nk

    klen[i1:i2], kdat[i1:i2][:] = self.set_kvec([ki1, ki2], [kf1, kf2], nk)

    ir += 1
    ik += nk

# G-K
    ki1 = 0.0
    ki2 = 0.0
    kf1 = 4*pi/3
    kf2 = 2*pi/3
#  kf1 = pi
#  kf2 = pi/2

    i1 = ik
    i2 = ik + nk

    klen[i1:i2], kdat[i1:i2][:] = self.set_kvec([ki1, ki2], [kf1, kf2], nk)

    return klen, kdat, nr, lett


# ********************************************************************
# ********************************************************************
#
# class BS_Tri
#
# ********************************************************************
# ********************************************************************
class BS_Tri_L3( Band_Structure ):
  """
  Simple triangular lattice
  """

# ********************************************************************
#
# __init__
#
# ********************************************************************
  def __init__( self ):
    """
    """
    self.nb = 1
    self.no = 3

    self.thop, self.eonsite = Init_Hops_NaxCoO2_Frank_L3()

#
# BZ basis vectors
#
    self.BZ_bas = 2 * pi * array( [[1, -1/sqrt(3)], \
                                   [0,  2/sqrt(3)]] )

#
# Super BZ (smaller BZ of the supercell)
#
    self.SBZ_bas = 2 * pi * array( [[1, -1/sqrt(3)], \
                                    [0,  2/sqrt(3)]] )
# 
#  Cluster basis vectors 
#  (defined in terms of a1, a2 - translation vectors of the original lattice)
# 
    self.sbas = array([[0.0, 0.0]])

    Band_Structure.__init__( self )


# ********************************************************************
#
# Hamiltonian
#
# ********************************************************************
  def Hamiltonian( self, kp ):
    """
    Evaluates the Hamiltonian for a given point 'k'
    """

    k1 = kp[0]
    k2 = kp[1]

    thop = self.thop
    eonsite = self.eonsite

    emat0 =  diag( eonsite, 0 )

    emat1 = -( thop[:,:,1] * exp( -1.0j * k1 ) + \
               thop[:,:,0] * exp( -1.0j * k2 ) + \
               thop[:,:,2] * exp( -1.0j * (k2 - k1) ) )
    emat1 += emat1.conj().T

    emat2 = -( thop[:,:,3] * exp( -1.0j * (k1 + k2) ) + \
               thop[:,:,4] * exp( -1.0j * (2*k2 - k1) ) + \
               thop[:,:,5] * exp( -1.0j * (2*k1 - k2) ) )
    emat2 += emat2.conj().T

    emat3 = -( thop[:,:,7] * exp( -2.0j * k1 ) + \
               thop[:,:,6] * exp( -2.0j * k2 ) + \
               thop[:,:,8] * exp( -2.0j * (k2 - k1) ) )
    emat3 += emat3.conj().T

    emat_full = emat0 + emat1 + emat2 + emat3
 
    return emat_full

# ********************************************************************
#
# High_Symmetry_kPoints
#
# ********************************************************************
  def High_Symmetry_kPoints(self, nk):
    """
    Generates lines along high-symmetry directions in k-space.
    """

    nr = 2

    klen = zeros(nk * nr, float64)
    kdat = zeros((nk * nr, 2), float64)

#
# Define 2 k-ranges
#
    ik = 0
    ir = 0
#
# Symbols for special k-points
#
#  lett = ['{/Symbol G}', 'K', 'M', '{/Symbol G}']
    lett = ['M', '{/Symbol G}', 'K']

# M-G
    ki1 = pi
    ki2 = pi
    kf1 = 0.0
    kf2 = 0.0

    i1 = ik
    i2 = ik + nk

    klen[i1:i2], kdat[i1:i2][:] = self.set_kvec([ki1, ki2], [kf1, kf2], nk)

    ir += 1
    ik += nk

# G-K
    ki1 = 0.0
    ki2 = 0.0
    kf1 = 4*pi/3
    kf2 = 2*pi/3
#  kf1 = pi
#  kf2 = pi/2

    i1 = ik
    i2 = ik + nk

    klen[i1:i2], kdat[i1:i2][:] = self.set_kvec([ki1, ki2], [kf1, kf2], nk)

    return klen, kdat, nr, lett


# ********************************************************************
# ********************************************************************
#
# class BS_Square
#
# ********************************************************************
# ********************************************************************
class BS_Square( Band_Structure ):
  """
  Simple square lattice
  """

# ********************************************************************
#
# __init__
#
# ********************************************************************
  def __init__( self ):
    """
    """
    self.nb = 1
    self.no = 1

    self.thop = 1.0

#
# BZ basis vectors
#
    self.BZ_bas = 2 * pi * array( [[1, 0], \
                                   [0, 1]] )

#
# Super BZ (smaller BZ of the supercell)
#
    self.SBZ_bas = 2 * pi * array( [[1, 0], \
                                    [0, 1]] )
# 
#  Cluster basis vectors 
#  (defined in terms of a1, a2 - translation vectors of the original lattice)
# 
    self.sbas = array([[0.0, 0.0]])

    self.svec = zeros( (2, 1), dtype=float )

    Band_Structure.__init__( self )


# ********************************************************************
#
# Hamiltonian
#
# ********************************************************************
  def Hamiltonian( self, kp ):
    """
    Evaluates the Hamiltonian for a given point 'k'
    """

    k1 = kp[0]
    k2 = kp[1]

    thop = self.thop

    emat_full = array([[-2.0 * thop * (cos(k1) + cos(k2))]])
 
    return emat_full

# ********************************************************************
# ********************************************************************
#
# class BS_Square_L2
#
# ********************************************************************
# ********************************************************************
class BS_Square_L2( Band_Structure ):
  """
  Simple square lattice, 2 orbitals
  """

# ********************************************************************
#
# __init__
#
# ********************************************************************
  def __init__( self ):
    """
    """
    self.nb = 1
    self.no = 2

    self.thop = 1.0

#
# BZ basis vectors
#
    self.BZ_bas = 2 * pi * array( [[1, 0], \
                                   [0, 1]] )

#
# Super BZ (smaller BZ of the supercell)
#
    self.SBZ_bas = 2 * pi * array( [[1, 0], \
                                    [0, 1]] )
# 
#  Cluster basis vectors 
#  (defined in terms of a1, a2 - translation vectors of the original lattice)
# 
    self.sbas = array([[0.0, 0.0]])

    Band_Structure.__init__( self )

# ********************************************************************
#
# Init_Hamiltonian
#
# ********************************************************************
  def Init_Hamiltonian( self, nk ):
    """
    Generate the k-mesh and Hamiltonian.
    """

    Band_Structure.Init_Hamiltonian( self, nk )

    self.svec = vstack( (self.sbas[0,:], self.sbas[0,:]) )

    print " svec:", self.svec.shape
    print self.svec

# ********************************************************************
#
# Hamiltonian
#
# ********************************************************************
  def Hamiltonian( self, kp ):
    """
    Evaluates the Hamiltonian for a given point 'k'
    """

    k1 = kp[0]
    k2 = kp[1]

    thop = array( [ [0.2, 0.05], [0.05, 0.2] ] )

    eonsite = array( [-0.0, 0.0] )

    emat_full = -2.0 * thop * (cos(k1) + cos(k2)) + diag( eonsite, 0 )
 
    return emat_full


# ********************************************************************
# ********************************************************************
#
# class BS_Kag_L3
#
# ********************************************************************
# ********************************************************************
class BS_Kag_L3( Band_Structure ):
  """
  Kagome lattice; 3-site basis.
  3-oribtals: a1g, e'g, e'g.
  """

# ********************************************************************
#
# __init__
#
# ********************************************************************
  def __init__( self ):
    """
    """
    self.nb = 3
    self.no = 3

#    self.thop, self.eonsite = Init_Hops_L3_deb()
    self.thop, self.eonsite = Init_Hops_NaxCoO2_Frank_L3()
#    self.thop, self.eonsite = Init_Hops_L3_v2_x33()

#
# BZ basis vectors
#
    self.BZ_bas = 2 * pi * array( [[1, -1/sqrt(3)], \
                                   [0,  2/sqrt(3)]] )

#
# Super BZ (smaller BZ of the supercell)
#
    self.SBZ_bas = pi * array( [[1, -1/sqrt(3)], \
                                [0,  2/sqrt(3)]] )
# 
#  Cluster basis vectors 
#  (defined in terms of a1, a2 - translation vectors of the original lattice)
# 
    self.sbas = array([[0.0, 0.0],
                       [1.0, 0.0],
                       [0.0, 1.0] ])

    Band_Structure.__init__( self )

# ********************************************************************
#
# Hamiltonian
#
# ********************************************************************
  def Hamiltonian( self, kp ):
    """
    Evaluates the Hamiltonian for a given point 'k'
    """

    k1 = kp[0]
    k2 = kp[1]

    thop = self.thop
    eonsite = self.eonsite

    emat0 = zeros((3, 3, 3, 3), complex64)

    emat0[0,0,:] =  diag( eonsite, 0 )
    emat0[1,1,:] =  emat0[0,0,:]
    emat0[2,2,:] =  emat0[0,0,:]

    emat1 = zeros((3, 3, 3, 3), complex64)
    emat2 = zeros((3, 3, 3, 3), complex64)
    emat3 = zeros((3, 3, 3, 3), complex64)

    emat1[1,0,:] = -( thop[:,:,1].T + thop[:,:,1] * exp( -2.0j * k1 ) )
    emat1[2,0,:] = -( thop[:,:,0].T + thop[:,:,0] * exp( -2.0j * k2 ) )
    emat1[2,1,:] = -( thop[:,:,2] + thop[:,:,2].T * exp( -2.0j * (k2-k1) ) )

    emat2[1,0,:] = -( thop[:,:,4] * exp( -2.0j * k2 ) + thop[:,:,4].T * exp( -2.0j * (k1-k2) ) )
    emat2[2,0,:] = -( thop[:,:,5] * exp( -2.0j * k1 ) + thop[:,:,5].T * exp( -2.0j * (k2-k1) ) )
    emat2[2,1,:] = -( thop[:,:,3] * exp( -2.0j * k2 ) + thop[:,:,3].T * exp( 2.0j * k1 ) )

    e3 = - ( thop[:,:,7] * exp( -2.0j * k1 ) + \
             thop[:,:,6] * exp( -2.0j * k2 ) + \
             thop[:,:,8] * exp( -2.0j * (k1 - k2) ) )

    emat3[0,0,:] = e3 + e3.conj().T
    emat3[1,1,:] = emat3[0,0,:]
    emat3[2,2,:] = emat3[0,0,:]

    emat0b = bmat( [[emat0[0,0,:], emat0[0,1,:], emat0[0,2,:]], \
                    [emat0[1,0,:], emat0[1,1,:], emat0[1,2,:]], \
                    [emat0[2,0,:], emat0[2,1,:], emat0[2,2,:]] ] )

    emat1b = bmat( [[emat1[0,0,:], emat1[0,1,:], emat1[0,2,:]], \
                    [emat1[1,0,:], emat1[1,1,:], emat1[1,2,:]], \
                    [emat1[2,0,:], emat1[2,1,:], emat1[2,2,:]] ] )

    emat2b = bmat( [[emat2[0,0,:], emat2[0,1,:], emat2[0,2,:]], \
                    [emat2[1,0,:], emat2[1,1,:], emat2[1,2,:]], \
                    [emat2[2,0,:], emat2[2,1,:], emat2[2,2,:]] ] )

    emat3b = bmat( [[emat3[0,0,:], emat3[0,1,:], emat3[0,2,:]], \
                    [emat3[1,0,:], emat3[1,1,:], emat3[1,2,:]], \
                    [emat3[2,0,:], emat3[2,1,:], emat3[2,2,:]] ] )
#
# Sort block by site indices rather than orbitals
#
# Permutations: 2 <-> 4, 3 <-> 7, 6 <-> 8
    Pm = eye(9)
# 2 <-> 4
    Pm[1,1] = Pm[3,3] = 0
    Pm[1,3] = Pm[3,1] = 1
# 3 <-> 7
    Pm[2,2] = Pm[6,6] = 0
    Pm[2,6] = Pm[6,2] = 1
# 6 <-> 8
    Pm[5,5] = Pm[7,7] = 0
    Pm[5,7] = Pm[7,5] = 1


    emat1b += emat1b.conj().T
    emat2b += emat2b.conj().T

    emat_full = dot( dot(Pm, emat0b + emat1b + emat2b + emat3b), Pm )
 
    return emat_full

# ********************************************************************
#
# Spec_Fun
#
# ********************************************************************
  def Spec_Fun( self, Chem_pot, k, en, Sigma_k_fun = None, unfold = True, 
                orbs = None ):
    """
    Evaluates the spectral function for a given k-point 'k'
    and energy mesh 'en'.

       * unfold = True : apply the up-folding procedure
         unfold = False: do not apply up-folding

    """

    return Band_Structure.Spec_Fun(self, Chem_pot, k, en, Sigma_k_fun, 
                                   unfold, orbs, site_block=False)

# ********************************************************************
#
# High_Symmetry_kPoints
#
# ********************************************************************
  def High_Symmetry_kPoints(self, nk):
    """
    Generates lines along high-symmetry directions in k-space.
    """

    nr = 2

    klen = zeros(nk * nr, float64)
    kdat = zeros((nk * nr, 2), float64)

#
# Define 2 k-ranges
#
    ik = 0
    ir = 0
#
# Symbols for special k-points
#
#  lett = ['{/Symbol G}', 'K', 'M', '{/Symbol G}']
    lett = ['M', '{/Symbol G}', 'K']

# M-G
    ki1 = pi
    ki2 = pi
    kf1 = 0.0
    kf2 = 0.0

    i1 = ik
    i2 = ik + nk

    klen[i1:i2], kdat[i1:i2][:] = self.set_kvec([ki1, ki2], [kf1, kf2], nk)

    ir += 1
    ik += nk

# G-K
    ki1 = 0.0
    ki2 = 0.0
    kf1 = 4*pi/3
    kf2 = 2*pi/3
#  kf1 = pi
#  kf2 = pi/2

    i1 = ik
    i2 = ik + nk

    klen[i1:i2], kdat[i1:i2][:] = self.set_kvec([ki1, ki2], [kf1, kf2], nk)

    return klen, kdat, nr, lett

# ********************************************************************
# ********************************************************************
#
# class BS_Kag
#
# ********************************************************************
# ********************************************************************
class BS_Kag( Band_Structure ):
  """
  Kagome lattice; 3-site basis
  """

# ********************************************************************
#
# __init__
#
# ********************************************************************
  def __init__( self ):
    """
    """
    self.nb = 3
    self.no = 1

# Frank's parametrization
    self.thop = -0.134
    self.thop2 = 0.028
    self.thop3 = 0.024
    self.thopz = 0.015
# Modified Frank's parametrization
#    self.thop = -0.152
#    self.thop2 = 0.020
#    self.thop3 = 0.014
#    self.thopz = 0.009

#
# BZ basis vectors
#
    self.BZ_bas = 2 * pi * array( [[1, -1/sqrt(3)], \
                                   [0,  2/sqrt(3)]] )

#
# Super BZ (smaller BZ of the supercell)
#
    self.SBZ_bas = pi * array( [[1, -1/sqrt(3)], \
                                [0,  2/sqrt(3)]] )
# 
#  Cluster basis vectors 
#  (defined in terms of a1, a2 - translation vectors of the original lattice)
# 
    self.sbas = array([[0.0, 0.0],
                       [1.0, 0.0],
                       [0.0, 1.0] ])

    Band_Structure.__init__( self )

# ********************************************************************
#
# Hamiltonian
#
# ********************************************************************
  def Hamiltonian( self, kp ):
    """
    Evaluates the Hamiltonian for a given point 'k'
    """

    k1 = kp[0]
    k2 = kp[1]

    emat = zeros( (3, 3), complex64 )
    ematp = zeros( (3, 3), complex64 )
    emat3 = zeros( (3, 3), complex64 )

    emat[0,1] = -(1 + exp( 2.0j * k1) )
    emat[1,0] = -(1 + exp(-2.0j * k1) )
    emat[0,2] = -(1 + exp( 2.0j * k2) )
    emat[2,0] = -(1 + exp(-2.0j * k2) )
    emat[1,2] = -(1 + exp( 2.0j * (k2 - k1)) )
    emat[2,1] = -(1 + exp(-2.0j * (k2 - k1)) )
    
    emat *= self.thop

    ematp[1,0] = -( exp( 2.0j * (k2 - k1)) + exp(-2.0j * k2) )
    ematp[2,0] = -( exp(-2.0j * (k2 - k1)) + exp(-2.0j * k1) )
    ematp[2,1] = -( exp( 2.0j * k1) + exp(-2.0j * k2) )
    ematp[0,1] = conjugate( ematp[1,0] )
    ematp[0,2] = conjugate( ematp[2,0] )
    ematp[1,2] = conjugate( ematp[2,1] )

    ematp *= self.thop2

    e3 = - 2 * self.thop3 * ( cos( 2*k1 ) + cos( 2*k2 ) + cos( 2*(k2 - k1) ) )

    emat3[0,0] = e3
    emat3[1,1] = e3
    emat3[2,2] = e3

    emat_full = emat + ematp + emat3

    return emat_full

# ********************************************************************
#
# High_Symmetry_kPoints
#
# ********************************************************************
  def High_Symmetry_kPoints(self, nk):
    """
    Generates lines along high-symmetry directions in k-space.
    """

    nr = 2

    klen = zeros(nk * nr, float64)
    kdat = zeros((nk * nr, 2), float64)

#
# Define 2 k-ranges
#
    ik = 0
    ir = 0
#
# Symbols for special k-points
#
#  lett = ['{/Symbol G}', 'K', 'M', '{/Symbol G}']
    lett = ['M', '{/Symbol G}', 'K']

# M-G
    ki1 = pi
    ki2 = pi
    kf1 = 0.0
    kf2 = 0.0

    i1 = ik
    i2 = ik + nk

    klen[i1:i2], kdat[i1:i2][:] = self.set_kvec([ki1, ki2], [kf1, kf2], nk)

    ir += 1
    ik += nk

# G-K
    ki1 = 0.0
    ki2 = 0.0
    kf1 = 4*pi/3
    kf2 = 2*pi/3
#  kf1 = pi
#  kf2 = pi/2

    i1 = ik
    i2 = ik + nk

    klen[i1:i2], kdat[i1:i2][:] = self.set_kvec([ki1, ki2], [kf1, kf2], nk)

    return klen, kdat, nr, lett



